function [t_Ws_L,S_Ws_L,...
          t_Ws_R,S_Ws_R,...
          t_Wu_R,S_Wu_R,...
          t_Wu_L,S_Wu_L] = PropagateManifold_Colinear_rk4(~,varargin)

%%

numvarargs = length(varargin);
% Number of optional arguments

if mod(numvarargs, 2) == 1
    error('Error: Please check name-value pair arguments');
end

name_arguments = varargin(1:2:end);
value_arguments = varargin(2:2:end);
% Initialize name-value arguments

for ii = 1:length(name_arguments)
    switch lower(name_arguments{ii})
        case 'poinitialstate'
            Xo = value_arguments{ii};
        case 'poperiod'
            tau = value_arguments{ii};
        case 'massratio'
            u = value_arguments{ii};
        case 'manifoldtrajectoryamount'
            numpoints = value_arguments{ii};
        case 'integrationtime'
            integrationTime = value_arguments{ii};
        case 'system'
            if value_arguments{ii} == 'SE'
                epsilon=5*10^(-6);
            elseif value_arguments{ii} == 'EM'
                epsilon=5*10^(-3);
            end
        case 'dt'
            dt = value_arguments{ii}
    end
end


%% Initial Guess
% northern

% tau = [1.39253390838970];
%[] Orbital period initial guess.
T = tau;
%[] Single orbital period

to = linspace(0,T,numpoints);
%[]

options=odeset('RelTol',1e-15,'AbsTol',1e-21);
%[] Increase the Ode 45 tol

[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S';
t = t';
%[] Integrate states.

Monodromy = State2STM(t,S,u);
%[] Monodromy matrix

[eig_vec,eig_val] = eigs(Monodromy);
%[] Find eigen values and vectors

eig_val = diag(eig_val);
%[] Get just the value.

[~,ind_stable] = min(eig_val);
[~,ind_unstable] = max(eig_val);
%[] Determine stable and unstable index

nu_stable = eig_vec(:,ind_stable);
nu_unstable = eig_vec(:,ind_unstable);
%[] Stable and unstable direction perturbation

InitialConditions_stable_outside = zeros(6,length(t));
InitialConditions_stable_inside = zeros(6,length(t));
InitialConditions_unstable_outside = zeros(6,length(t));
InitialConditions_unstable_inside = zeros(6,length(t));

%[] Pre allocate initial conditions matrix

% Loop for creating initial conditions for manifold trajectories.
parfor ii = 1:(length(t)-1)
    
    t_STM = [t(1),t(ii+1)];
    %[] Determine the state transition matrix at desired time T
    
    Pert_unstable = State2STM(t_STM,Xo,u)*nu_unstable;
    Pert_stable = State2STM(t_STM,Xo,u)*nu_stable;
    %[] calculate the perturbance needed for unstable and stable manifold
    %trajectories
    
    Norm_Pert_unstable = Pert_unstable / norm(Pert_unstable);
    Norm_Pert_stable = Pert_stable / norm(Pert_stable);
    %[] Normalize the perturbing force.
    
    if Norm_Pert_stable(1)<0 % If x perturbation direction is negative
        
        InitialConditions_stable_outside(:,ii)      = S(:,ii) + epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) - epsilon * Norm_Pert_stable;
%         disp('a')
    else
        InitialConditions_stable_outside(:,ii)      = S(:,ii) - epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) + epsilon * Norm_Pert_stable;
%         disp('b')
    end
    
    if Norm_Pert_unstable(1)<0 % If x perturbation direction is negative
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) + epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) - epsilon * Norm_Pert_unstable;
%         disp('a')
    else
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) - epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) + epsilon * Norm_Pert_unstable;
%         disp('b')
    end
    
    %[] Determine the initial states for stable manifold's initial states.
    
end
T_manifold_stable_outside = integrationTime:-dt:0; % Integrating backward in time
%[] Integration time for stable manifold coming from outside

T_manifold_unstable_inside = 0:dt:integrationTime;% Integrating forward in time
% options_SecondBodyCrossing=odeset('RelTol',1e-6,'AbsTol',1e-9,'Events',@(t,S)HaltProp(t,S,u));

N=length(T_manifold_unstable_inside)-1;




for ii = 1:length(InitialConditions_stable_outside)-1
    S_Ws_L{ii}(:,1) = InitialConditions_stable_outside(:,ii);
    S_Ws_R{ii}(:,1) = InitialConditions_stable_inside(:,ii);
    S_Wu_R{ii}(:,1) = InitialConditions_unstable_inside(:,ii);
    S_Wu_L{ii}(:,1) = InitialConditions_unstable_outside(:,ii);

    for jj = 1:N
        S_Ws_L{ii}(:,jj+1) = f_rk4(@dynamics,0,S_Ws_L{ii}(:,jj),[u],-dt);
        S_Ws_R{ii}(:,jj+1) = f_rk4(@dynamics,0,S_Ws_R{ii}(:,jj),[u],-dt);
        S_Wu_R{ii}(:,jj+1) = f_rk4(@dynamics,0,S_Wu_R{ii}(:,jj),[u],dt);
        S_Wu_L{ii}(:,jj+1) = f_rk4(@dynamics,0,S_Wu_L{ii}(:,jj),[u],dt);
    end
    t_Ws_L{ii} = T_manifold_stable_outside;
    t_Ws_R{ii} = T_manifold_stable_outside;
    t_Wu_R{ii} = T_manifold_unstable_inside;
    t_Wu_L{ii} = T_manifold_unstable_inside;
end






end

function dS  = dynamics(t,S,u)
    x = S(1);
    y = S(2);
    z = S(3);
    vx = S(4);
    vy = S(5);
    vz = S(6);
    %[] saperate states into each components

    r1 = ((x+u)^2+y^2+z^2)^(1/2);
    r2 = ((x+u-1)^2+y^2+z^2)^(1/2);
    %[] Range of satellite wrt first and second primary bodies.
    
    n = length(S);
    %[] Length of the state vector!
    
    dS = zeros(n,1);
    %[] Pre allocate memory

    dS(1:3) = S(4:6);
    %[] Derivative of the position are the velosity at that time.
    
    dS(4) = x  -  (1-u)*(x+u)/r1^3  -  u*(x+u-1)/r2^3 + 2*vy;
    dS(5) = y  -  (1-u)*y/r1^3 - u*y/r2^3 - 2*vx;
    dS(6) = -(1-u)*z/r1^3-u*z/r2^3;
    %[] Acceleration in the rotating CRTBP frame.

end

function x_i = f_rk4(fun, t, x, par, dt)



k1 = fun(t, x, par)*dt;
k2 = fun(t + 0.5*dt, x + 0.5*k1, par)*dt;
k3 = fun(t + 0.5*dt, x + 0.5*k2, par)*dt;
k4 = fun(t + dt, x + k3, par)*dt;

x_i = x + 1/6*(k1 + 2*k2 + 2*k3 + k4);
end